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ABSTRACT 

This  thesis  deals  with  the  enhancement  of  video  images 
degraded  by  turbid  water  viewing  conditions.  An  algorithm  by  Peli 
and  Lim  has  been  used  with  some  success  for  enhancement,  but  it 
was  found  to  accentuate  noise.  The  thesis  examines  a  combination 
of  the  Peli  and  Lim  algorithm  with  three  approaches  to 
enhancement. 

First,  a  Short  Space  Spectral  Subtraction  algorithm  which 
performs  the  restoration  in  the  density  domain,  using  an  estimate 
for  the  power  spectrum  of  the  given  data  set.  The  degraded  image 
is  divided  into  many  subimages  and  each  subimage  is  restored 
separately  and  then  combined. 

Next,  an  algorithm  for  Image  Enhancement  and  Noise  Filtering 
by  Use  of  Local  Statistics,  which  uses  the  assumption  that  the 
sample  mean  and  variance  of  a  pixel  is  equal  to  the  local  mean 
and  variance  of  all  pixels  within  a  fixed  range  surrounding  it, 

Finally,  a  median  filter  for  noise  reduction  ,  where  a  given 
pixel  of  a  degraded  image  is  replaced  by  the  median  of  the  pixel 
values  in  a  window  surrounding  it. 

Combination  of  the  algorithms  are  applied  to  degraded  images, 
and  the  results  are  compared  and  discussed,  in  each  case.  It  was 
found    that    noise    smoothing    can    be    achieved    with    the    spectral 


subtraction    algorithm,     and    that    the    local    statistics     technique 
yielded  very  good  contrast  enhancement. 
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I.   INTRODUCTION 

In  this  thesis,  some  algorithms  for  the  enhancement  of  video 
images  degraded  by  turbid  water  viewing  conditions  are 
implemented.  The  images  used  to  illustrate  the  processes,  were 
recorded  during  the  recovery  procedure  in  a  torpedo  testing  range. 

The  recovery  operation  is  carried  on  undersea  with  the  aid  of 
digging  equipment  and  underwater  video  cameras  equipped  with 
strong  artificial  lights.  The  equipment  is  controlled  and  the 
operation  is  monitored  on  video  monitors  at  the  surface  onboard 
the  recovery  vessel.  The  recovery  equipment  in  its  attempt  to  dig 
out  the  torpedo,  stirs  up  sediment  wich  visually  obscures  the 
object  of  interest  and  impedes  the  operation.  To  date  no  special 
techniques  have  been  used  operationally  to  process  or  enhance  the 
video  image  before  display. 

When  it  is  of  interest  to  "improve"  an  image,  there  are  two 
broad  types  of  image  manipulation  processes  that  cover  all 
operations  performed  to  get  such  improvement:  image  restoration 
and  image  enhancement. 

The  goal  of  image  enhancement  is  to  process  a  degraded  image 
so  that  the  result  is  more  suitable  than  the  original  image  for  a 
specific  application  or  aids  the  human  analyst  in  the  extraction 
and  interpretation  of  pictorial  information  [Ref.  l] .  Image 
restoration,    on   the   other   hand,    is   intended   to   bring   a   degraded 


image  back  to  an  ideal  degradation-free  image  as  closely  as  possible 
[Ref.  2].  This  makes  the  image  enhancement  problem  a 
subjective  one,  the  improvement  of  the  image  appearance  to  the 
human  viewer  is  highly  dependent  on  the  viewer  himself  and  on 
the  application.  What  is  "good"  for  one  person  or  application  is  not 
necessarily  "good"  for  another  one.  Consequently,  the  suitability  of 
the  processed  image  for  a  specific  application  makes  image 
enhancement  techniques  very  much  problem-oriented. 

The  degradation  process  can  result  from  causes  such  as,  an 
imperfect  photographic  process,  imperfect  display  devices,  poor 
contrast  due  to  environmental  conditions,  different  forms  of  noise 
(channel,  quantization,  salt  and  pepper)  and  others  [Ref.  3].  In 
the  same  manner,  techniques  for  enhancement  are  dependent  of 
the  original  degradation  process.  Processes  yielding  satisfactory 
results  for  one  type  of  degradation  are  not  necessarily  suitable  for 
another. 

The  interest  of  this  thesis  is  on  examining  different  techniques 
to  enhance  images  that  are  degraded  by  environmental  conditions 
that  create  a  great  lack  of  contrast  in  the  recorded  images.  This 
work  is  a  continuation  of  previous  research  begun  on  the  subject, 
which  included  the  implementation  of  an  adaptive  filtering  [Ref.  4] 
algorithm  for  contrast  enhancement.  This  algorithm  was  applied  to 
a  set  of  underwater  images  recorded  in  turbid  water  viewing 
conditions.  The  images  are  low  contrast  degraded  by  noisy 
background.     As    will    be    explained     in     the     next     chapter,     the 
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algorithm  yielded  very  good  contrast  enhancement  but,  seemed  to 
accentuate  the  noise.  The  purpose  of  this  work  is  to  investigate 
some  noise  smoothing  techniques  to  be  applied  to  the 
contrast-enhanced  images,  an  alternative  contrast  manipulation 
scheme,  and  to  compare  the  results  obtained  to  those  found  with 
the  algorithm  already  implemented.  Besides  the  images  included  in 
the  experimental  results  of  this  thesis,  many  more  images 
resulting  from  other  combinations  and  variations  of  the  processes 
and  different  settings  of  the  parameters  were  obtained.  We  have 
included  here  only  those  thought  to  be  most  representative  of  the 
results  obtained. 


II.     OVERVIEW  OF  PREVIOUS  WORK 

A.     DESCRIPTION  OF  THE  ADAPTIVE  FILTERING  ALGORITHM 

An  adaptive  filtering  algorithm  for  image  enhancement  by 
Peli  and  Lim  [Ref.  4]  has  been  used  for  contrast  enhancement  of 
images.  The  images  are  typifyed  by  a  large  dynamic  range,  and 
recorded  over  a  medium,  characterized  by  a  much  smaller 
dynamic  range,  which  causes  those  regions  of  the  images  with 
very  high  or  very  low  luminance  to  be  poorly  represented.  This 
algorithm  modifies  the  local  lunninance  mean  of  an  linage  and 
controls  the  local  contrast  as  a  function  of  the  local  luminance 
mean  of  the  image. 


f(nl,n2) 
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Figure  2.1.     Adaptive  Filtering  for  Image  Enhancement 

In    Figure    2.1    a    block    diagram    of    the    Adaptive    Filtering 
algorithm     is     shown  In     the     figure,      f(ni,n2)     denotes     the 
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unprocessed  digital  image,  and  fL(ni,n2), which  denotes  the  local 
luminance  mean  of  f(ni,n2),  is  obtained  by  low- pass  filtering 
f(ni,n2).  The  sequence  fH(ni,n2)  which  denotes  the  local  contrast 
is  obtained  by  subtracting  fL(ni,n2)  from  f(ni,n2).  The  local 
contrast  is  modified  by  multiplying  fH(ni,n2)  with  k(fL),  a  scalar 
which  IS  a  function  of  fL(ni,n2).  The  modified  contrast  is  denoted 
by  f'H(ni,n2).  The  specific  functional  form  of  k(fL)  depends  on  the 
particular  application  under  consideration,  a  value  of  k(fL)  >  1 
represents  local  contrast  increase  while  k(fL)  <  1  represents  local 
contrast  decrease  The  local  lunninance  mean  is  modified  by  a  point 
non-linearity  and  the  modified  local  luminance  mean  is  denoted  by 
f'L(ni,n2).     The     specific     non-linearity    chosen     depends     on     the 

particular  application  under  consideration,  and  in  most  application 
problems  the  non-linearity  is  chosen  so  that  the  overall  dynamic 
range  of  the  resulting  image  is  approximately  the  same  as  the 
dynamic  range  of  the  recording  medium.  The  modified  local 
contrast  and  local  luminance  mean  are  then  combined  to  obtain 
g(ni,n2),   the  processed  image. 

B.     DISCUSSION  OF  RESULTS 

In  this  section  the  application  of  the  algorithm  described  in  the 
previous  section  to  the  images  shown  in  Figure  2.2  is  presented . 
Figure  2.2  shows  two  low-contrast  noise-degraded  images  of  512  X 
512  pixels  with  each  pixel  represented  by  8  bits.  Note  that  the 
small    dynamic  range    of    the    lunninance    of    both   images    make 


11 


(b) 

Figure  2.2    Original  Degraded  Images. 
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difficult  the  appreciation  of  subtle  details  in  thenn.  In  this  case,  it 
is  desired  to  increase  the  local  contrast  in  the  images  to  bring  back 
the  lost  details.  To  achieve  this,  k(fL)  is  chosen  as  shown  in  Figure 
2.3  (a)  and  (b)  for  the  images  of  Figure  2.2(a)  and  (b), 
respectively.  For  the  case  of  Figure  2.2(a),  the  function  is  chosen 
to  increase  the  local  contrast  for  the  low  luminance  regions, 
maintaining  it  relatively  constant  for  medium  ranges  of  local 
luminance,  and  again  enhance  it,  for  the  brighter  regions  of  the 
degraded  image.  The  results  are  shown  in  Figure  2.4(a),  where  a 
relative  contrast  increase  can  be  observed.  This  is  especially 
evident  in  the  "fish"  shape  close  to  the  lower  right  corner  of  the 
image.  Also  note  the  increase  in  background  noise  resulting  in  the 
textured  background  observed  in  the  contrast-enhanced  image. 
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Figure  2.3.     K(fL)  Function  for  Processing  Images  in  Figure  2.2 
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(b) 
Figure  2.4.   Images  in  Figure  2 . 2  Processed  with  Adaptive  Filtering 
for  Image  Enhancement  Algorithm. 
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In  the  case  of  the  image  of  Figure  2.2(b),  it  is  desired  to 
retain  the  original  local  contrast  in  the  darker  regions,  but 
enhance  it  in  the  regions  corresponding  to  medium  and  high  local 
luminance  means.  This  explains  the  particular  shape  of  the  k(fL) 
function  chosen  to  process  this  image.  In  the  results,  shown  in 
Figure  2. 4(b), some  details  difficult  to  be  distinguished  in  the 
un-processed  image,  are  now  more  evident.  Among  them  are  the 
particular  shape  of  the  object  near  the  center  of  the  image  and 
the  form  of  the  ocean  floor  below  the  left-most  object.  Again,  it 
can  be  observed  that  there  is  an  increase  in  the  background  noise 
throughout  the  contrast-enhanced  image. 
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III.   DESCRIPTION  OF  THE  ALGORITHMS 

In  this  chapter,  the  theorethical  foundations  and  procedures 
for  implementation  of  the  different  algorithms  used  m  the 
noise-filtering  and  contrast  enhancement  problem  are  presented. 
We  start  with  a  description  of  noise-snaoothing  techniques  to  be 
used  as  post-processing  after  the  adaptive  filtering  algorithm.  Then 
an  alternative  method  for  the  contrast  enhancement  is  presented. 
Finally,  the  use  of  median  filtering  for  noise-smoothing  is 
described 

A.     SHORT  SPACE  SPECTRAL  SUBTRACTION 

This  algorithm,  developed  by  J.  S.  Lim  [Ref  5]  performs  the 
restoration  in  the  frequency  domain,  using  a  converging  solution 
for  the  power  spectrum  of  a  given  data  set  that  does  not  depend 
on  the  original  estinnate.  The  procedure  is  based  on  two  basic 
assumptions  First,  each  part  of  an  image  generally  differs 
sufficiently  from  other  parts  so  that  the  image  cannot  be  modeled 
by  a  stationary  random  process.  Second,  the  power  spectrum  of 
the  restored  linage  is  estimated  by  the  spectral  subtraction  of  the 
additive  noise  spectra  from  the  degraded  image. 

1.      Spectral  Subtraction 

Frequency  domain  techniques  for  image  enhancement   use 
iterative     procedures     for     estimating     the     image     power     spectral 
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density,  continuing  until  certain  fidelity  criteria  are  met.  The 
density  so  obtained,  is  dependent  on  both  the  initial  estimate  for 
starting  the  iterative  procedure  and  the  stopping  criterion.  Some 
methods  converge  more  rapidly  than  others,  and  sometimes  first 
and  second  order  statistics  of  the  degrading  function  must  be 
known . 

Consider  a  model  of  image  degradation  as  presented  in 
Figure  3.1.  In  this  model  f(ni,n2)  represents  a  digital  image, 
b(ni,n2)  represents  a  linear  space-invariant  point  spread  function, 
and  d(ni,n2)  represents  an  additive  noise  component.  Thus, 
g(ni,n2)  represents  a  blurred  version  of  f(ni,n2)  and  y(ni,n2)  is 
the  degraded  image. 


f(nl,n2) 


g(nl,n2) 


y(nl,n2) 


b(nl,n2) 


d(nl,n2) 


Figure  3.1,     Model  of  Image  Degradation 

An  image  restoration  system  is  shown  in  Figure  3.2.  It 
can  be  seen  that  the  goal  of  a  restoration  system  is  to  process  a 
degraded  image  y(ni,n2)  through  a  noise  reduction  system  to 
estimate    g(ni,n2)    which    in    turn    is    processed    by    a    deblurring 
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system  to  get  f(ni,n2),   an  estimate  of  the  original  blur-free  noise- 
free  image. 


f(nl,n2) 

y(nl,n2^ 

NOISE  REDUCTION 
SYSTEM 

g(nl,n2)^ 

DEBLURRING 
SYSTEM 

w 

Figure  3.2.     Image  Restoration  System 

The  purpose  of  this  algorithm  is  to  provide  a  model  for 
the  Noise  Reduction  System,  as  presented  in  the  Figure.  A  general 
model  for  the  restoration  filter  is  given  by 


H((i) 


\'V  '  { 


P  (d)  ,0)  ) 
9^    1     2' 


P  (0)  ,(0  )  +  a  •  P^((i)  ,(i)  ) 


} 


P 


(3.1) 


where  Pg ((1)1,(1)2)  and  Pd(wi,<^)  represent  the  power  spectra  of 
g(ni,n2)  and  d(ni,n2)  respectively,  a  and  p  are  constants.  If  P  is 
unity,   H ((1)1,  (1)2)  corresponds  to  the  parametric  Wiener  filter,   and  if 

a  also  IS  unity,   it  reduces  to  the  standard  Wiener  filter  [Ref.   1,2]. 
Other  values  for  a  and  p  define  different  filtering  techniques,    such 

as  power  spectrum  filtering  and  geometrical  mean  filtering.  One 
commonly  used  approach  in  the  implementation  of  H  (0)1,002)   is  an 

iterative  procedure  which  begins  with  an  initial  estimate  of 
Pg((x>i,G02)   and  then  iteratively  estimates  g(ni,n2)   and  Pg(coi,co2) 
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until  a  converging  solution  or  a  desirable  performance  is  achieved. 
This  iterative  procedure  may  be  computationally  undesirable. 

The  Spatial  Subtraction  approach  analytically  obtains  a 
converging  solution,  by  estimating  Pg (001,002)  from  g(ni,n2)  as 
(i/k)«|G  (001,002)1^  and  using  a  value  for  ^  of  V2.  The  complete 
solution    for  this  problem,   is  given  by 


|g(w. 0)^)1   =  ^  |y(w.(o^)  1'- a.k.p^(w.w  )     (3.2) 


^  G(w.w  )  =  ^  Y(w.w  )  (3.3) 


for      |Y(W.W)p^  a.k.P^(CO^.W  ) 
and  zero  otherwise 

Here  G  ((1)1,0)2)  represents  the  discrete  space  Fourier  Transform  of 
g(ni,n2),  ¥(0)1,(1)2)  ^s  the  discrete  space  Fourier  Transform  of 
y(ni,n2),  Pd(<«ii)^)  represents  the  power  spectrum  of  d(ni,n2), 
the  additive  noise  component  of  the  degradation  model,  and  k  is  a 
scaling  factor  that  normalizes  the  power  and  energy  spectral 
densities.  The  phase  of  g(ni,n2)  is  estimated  by  the  phase  of 
y(ni,n2)  and  the  transform  magnitude  of  g(ni,n2)  is  estimated  by 
a  particular  form  of  spectral  subtraction. 
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2.     Short  Space  Implementation 

As  mentioned  before,  the  basic  assunnption  is  that  each 
part  of  an  image  f(ni,n2)   generally  differs  sufficiently  from  other 

parts  so  that  it  cannot  be  modelled  by  a  stationary  random  field. 
Thus  a  short  space  implementation  seems  appropiate.  The 
technique  consists  of  dividing  the  original  image  into  many 
subimages  and  restoring  each  subimage  separately.  The  procedure 
involves  the  application  of  a  short  space  window  function 
Wij(ni,n2)   to  overlapping  portions  of  the  degraded  image  y(ni,n2) 

so  that: 

y(n^,n2)«w^  .^(n^,!!^)   =   {g(n^,n2)+d(n^,n2)  1    •   ▼i^-jlnj/ng)    '^^^^ 
or,   equivalently. 

The  noise  reduction  system  is  then  applied  to  yij(ni,n2) 
to  recover  gij(ni,n2).  The  overall  full-size  restored  image  g(ni,n2) 
is  given  by 


i=0     j=0 
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The  window  function  is  desired  to  be  a  smooth  function  to 
avoid  some  possible  discontinuities  that  may  appear  at  the 
subimage  boundaries  in  the  processed  image  and  must  satisfy: 


2K      2C 


2  2  ^,3(^'^2)  =  ^ 

1-0    j-0 
for  all  Hj ,n2  of  interest 


(3.7) 


_/ 
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w^Cn^)   w^Cn^-L)    W2(n2-2L) 

/ 
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n 

▼  (n^,!!^)    =  w^(nJ.W2(n2) 

2 

Figure  3.3.   2-D  Separable  Triangular  Window 

The  condition  in  Equation  3 . 6  guarantees  that  an  image 
can  be  reconstructed  from  its  subimages.  To  ensure  the  smoothness 
of  the  window,  a  2-D  separable  triangular  window  similar  to  that 
shown  m  Figure  3 . 3  will  be  used .  Each  window  overlaps  its 
neighboring  window  by  half  the  window  length  in  each  dimension. 
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B.     CONTRAST  ENHANCEMENT  BY  USE  OF  LOCAL  STATISTICS 

This  algorithm  presents  a  computational  technique  for  contrast 
enhancement  on  a  two-dimensional  image  array  based  on  the  local 
mean  and  variance.  The  algorithm,  as  developed  by  J.  Lee  [Ref. 
6]  uses  the  basic  assumption  that  the  sample  mean  and  variance 
of  a  pixel  is  equal  to  the  local  mean  and  variance  of  all  pixels 
within  a  fixed  range  surrounding  it. 

Let  Xij  be  the  brightness  of  a  pixel  (i,  j)     in  a  two  dimensional 

NXN   image.    The   local   mean   and   variance   are   calculated   over   a 
(2n  +l)x(2m  +l)  window.   The  local  mean  is  defined  as 


i-fn 


■°'.r(2n.lW2m*l)  X    ,4-  ''-i  '^^' 


)(2m*l)  ^^ 


i-n  l"j-ffl 


and  the  local  variance  is  defined  as 


"'.r,2n.iw2M.i) ,?  X  '^.r^.j''         <^^> 


The  algorithm  is  designed  such  that  a  pixel  Xjj  will  maintain 
its  local  mean,  and  yet  permit  its  variance  to  be  modified  by  a 
constant  factor  times  its  original  variance.  The  procedure  is  defined 
by 
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where  k,  the  gain,  is  the  ratio  of  new  local  standard  deviation  to 
the  original  standard  deviation.  This  approach  has  the 
computational  advantage  that  the  local  variance  v^ j  is  not 
required  and  only  the  local  mean  mij  needs  to  be  computed.  If 
we  have  k  >  1,  the  image  will  be  sharpened  as  if  acted  upon  by  a 
high-pass  filter.  If  0  <  k  <  1,  the  image  will  be  smoothed  as  if 
passed  through  a  low-pass  filter.  In  the  extreme  case, one  has  k  = 
0  and  x'ij  is  equal  to  its  local  mean  mij  . 


X 


>.j 


LOCAL 
MEAN 


Figure  3.4.     Contrast  Enhancement  by  Use  of  Local  Statistics 

Figure  3.4  shows  a  block  diagrana  of  the  algorithm.  The  local 
mean  is  evaluated  over  a  (2n+l)  (2m+l)  window  and  subtracted 
from   the   original   brightness   of   the    pixel,    Xij.    The    difference    is 
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then  multiplied  by  the  standard  deviation  ratio  to  enhance  the 
local  contrast.  Finally,  this  result  is  added  to  the  local  mean 
yielding  the  reconstructed  pixel  value  x'ij. 

C.     NOISE  FILTERING  BY  USE  OF  LOCAL  STATISTICS 

This  algorithm,  also  developed  by  J.Lee  [Ref.  6]  is  actually  a 
variation  of  the  previous  algorithm  adapted  for  the  noise  filtering 
problem.  In  this  algorithm,  the  a  prion  mean  (or  variance)  of 
the  estimated  image  is  calculated  as  the  difference  between  the 
mean  (or  variance)  of  the  noise  corrupted  image  and  the  mean 
(or  variance)  of  the  noise  itself.   Let  Zij  be  the  degraded  pixel  Xjj. 

Then  the  degraded  pixel  can  be  modelled  as  the  sum  of  the 
noise-free  pixel  and  a  white  noise  sequence.   That  is. 


^,j  =  ^i,j^  Vj  (^^^) 


where  Wij  is  the  white  random  sequence  with  E[wij  ]=0.    Define 
the  a  prion    mean  and  variance  of  Xj j  to  be 


and 


Oij=  E[(x,^j-x,^^)']=  E[(z^  .,-z-^^^)'-0j  (3.13) 

respectively.   The  estimated  pixel  value  x^ j  is  computed  by 
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where  the  gain  is  given  by 


(3.14) 


i'J        1 


(3.15) 
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Figure  3.5.     Noise  Filtering  by  Use  of  Local  Statistics 

Figure   3 . 5   shows  a  block  diagram   for   the   implementation   of 
the  algorithnn.   Starting  with  the  degraded  pixel  z^j.the  local  mean 

is   evaluated   and    subtracted    to    obtain    the    local    contrast.    Then, 
both    Qij    and    k^ j    can    be    evaluated,    and    the    value    of    k    thus 

obtained  is  multiplied  by  the  local  contrast.    This  in   turn   is  added 
to  the  local  mean  to  get  the  restored  pixel  x^ j. 

One    innportant    paranneter    is    the    size    of    the    window    over 
which  the  local  mean  is  estimated.    If  the  window  is  too  small,    the 
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noise  filtering  algorithm  is  not  effective.   If  the  window  is  too  large, 
subtle  details  of  the  image  will  be  lost  in  the  filtering  process. 

D.     MEDIAN  FILTERING 

Median  Filtering  is  a  non-linear  signal  processing  technique 
originally  developed  by  J.Tukey  [Ref.  7],  that  is  useful  for  noise 
suppression  m  images.  In  one-dimensional  form  the  median  filter 
consists  of  a  sliding  window  encompassing  an  odd  number  of  pixels. 
The  center  pixel  in  the  window  is  replaced  by  the  median  of  the 
pixels  in  the  window.   The  median  of  a  discrete  sequence  ai,a2,  .    . , 

a^  for  M  odd  is  that  member  of  the  sequence  for  which  (M-l)/2 
elements  are  smaller  or  equal  in  value,  and  (M-l)/2  elements  are 
larger  or  equal  in  value  [Ref.  2].  For  example,  if  the  ordered 
values  of  the  pixels  within  a  window  are  80,  90,  200,  110,  120, 
the  center  pixel  would  be  replaced  by  the  value  110,  which  is  the 
middle  point  of  the  sorted  sequence  80,90,110,120,200. 

The  concept  of  median  filter  is  easily  extended  to  two 
dimensions  by  utilizing  a  two-dimensional  window  of  some  desired 
shape  such  as  a  rectangle  or  a  discrete  approximation  to  a  circle. 
A  two-dimensional  LXL  median  filter  will  provide  a  greater  degree 
of  noise  suppression  than  sequential  horizontal  and  vertical 
processing  with  LXl  median  filters,  but  two-dimensional  processing 
also  results  in  greater  signal  suppression. 

The  median  filter  has  been  found  to  be  more  effective  than  a 
linear    filter    for    smoothing    images   with    spiky    noise    degradations 
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because  of  extrema  rejection  by  the  median   [Ref,    8]  .Furthermore, 
the  median  filter  preserves  monotonic  step  edges,    that  is,    it  does 
not  blur  sharp  edges  as  a  linear  low-pass  filter  would. 

Another  interesting  aspect  of  median  filtering  that  has  been 
studied  [Ref.  9],  is  that  of  the  convergence  of  a  filtered  image  to 
what  has  been  called  a  root  signal  .  A  root  is  a  signal  which  is 
invariant  under  filtering  by  a  particular  median  filter.  The 
technique  consists  in  making  succesive  passes  of  a  noise-degraded 
image  thru  a  median  filter  until  an  image  corresponding  to  a  root 
signal  is  achieved.  Definition  of  a  root  signal  follows,  where  a 
window  width  of  2N+1  is  used. This  definition  uses  the  following 
ideas: 

1)  A  constant  neighborhood  is  a  region  of  at  least  N+1  consecutive 
identically  valued  samples. 

2)  An  edge  is  an  increasing  or  decreasing  sequence  of  samples 
which  is  immediately  preceded  and  followed  by  constant 
neighborhoods.  An  edge  cannot  contain  any  constant 
neighborhood. 

3)  An  impulse  is  a  sequence  of  at  most  N  consecutive  samples 
whose  values  are  different  from  those  of  the  two  surrounding 
regions,  the  two  surrounding  regions  are  identically  valued 
constant  neighborhoods. 

A    signal    is    a    root   if    and    only    if    it    contains    only    edges 

alternating   with    constant    neighborhoods    [Ref.     10].     The    median 

filtering  preserves  edges  and  constant  neighborhoods  but  eliminates 

impulses.   The  technique  is  based  on  the  idea  that  passing  a  noise 

corrupted  root  signal  through  the  median  filter  a  sufficient  number 
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of  times  will  produce  another  root  signal,  which  is  usually  "close" 
to  the  original  root  signal.  In  reference  [Ref.  lO],  the  concept  of 
"closeness"  and  a  measure  of  the  number  of  runs  to  be  used  are 
calculated  as  a  function  of  the  image  and  window  sizes. 

In  this  thesis,    a  median  filter  included  in   the   Spider   library 
[Ref.    11]   IS  used,    which   is   the   implementation   of    the   algorithm 
suggested    by     Huang,     Yang     and     Tang     [Ref.      12].      It     consists 
primarily  of  an  efficient  way  of  updating  the  histogram  each  time 
that  a  given  pixel  is  replaced  by  the  median  within  a  window     The 
algorithm  is  described  briefly  for  a  3x3  window: 
Step  1 :         Obtain  the  histogram  in  the  first  window  [Figure    3 . 6 
(a)]  and  find  the  median  MDN.   Next, count  the  pixels 
with  values  smaller  than  the  median  of  the  window. 
LTMDN  is  that  number. 
Step  2 :         Shift  the  window  right  by  one  pixel  [Figure  3 . 6  (b)  ] 
and  update  the  histogram  and  LTMDN.   First  decre- 
ment the  counted  values  of  the  histogram  equivalent 
to  gray  level  values  {gJ (a) ,  gJ fd) ,  gJ (gj  )  of  pixels  a,    d 
and  g.  That  is: 

mS'T{gl(a)  ]  =  m^l\gl(a)  ]-l 
msi\gl(d)  ]  =  msi{gl(d)  ]-l 
mzi{gl(g)  ]  =  ^\si{gl(g)  ]-l 
LTMDN  is  updated  as  follows: 

IF   gl(a)    <  MDN  ,   LTMDN  =  LTMDN- 1 
IF   g/{yj    <  MDN  ,   LTMDN  =  LTMDN- 1 
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IF    gl(g)    <  MDN  ,   LTMDN  =  LTMDN-1 
In  the  same  way  for  pixels  c',f',   and  i', 
msT[g/(c'J  ]  =  msT[g/{c'J  ]+l 
msT[g/{y'J  ]  =  msT[gJ(?'J  ]+l 
msT[gJ{j"J  ]  =  msT[g/{}'J  ]+l 
IF   g/{c'J    <  MDN  ,   LTMDN  =  LTMDN+1 
IF   gJ^/'J    <  MDN  ,   LTMDN  =  LTMDN+1 
IF   gJfj'J    <  MDN  ,   LTMDN  =  LTMDN+1 


Figure  3.6.   Window  for  Median  Filtering 

Step  3:         Update  MDN  in  the  previous  window,   and  obtain  a 
median  in  the  new  window.   Let  ITH=  (number  of 
pixels  in  the  window) /2 
-  If  LTMDN  >  ITH,   then 

LTMDN  =  LTMDN  -  HIST  [MDN] 
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MDN  =  MDN-1 
Repeat  this  until  LTMDN  <  ITH  is  obtained. 
-  If  LTMDN  <.  ITH,   then 

MDN  =  MDN+1 

LTMDN  =  LTMDN  +  HISt[MDN] 
Repeat  this  until  LTMDN  +  HISt[MDN+1]  >  ITH  is 
obtained.   This  MDN  is  the  median  for  the  current 
window  and  substitute  it  in  the  output  image. 
Step  4:         End  when  one  line  is  finished.  Go  to  Step  2. 
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IV.   EXPERIMENTAL  RESULTS 

In  this  chapter,  experimental  results  for  different  combinations 
of  the  algorithms  described  are  presented  and  discussed.  Basically, 
the  adaptive  filtering  algorithm  [Ref.  4]  was  used  for  contrast 
enhancement,  as  described  in  Chapter  II.  As  previously  stated, 
this  algorithm  yielded  very  good  contrast  enhancement,  but  also 
tended  to  accentuate  the  noise.  Thus,  methods  for  noise  reduction, 
such  as  short  space  spatial  subtraction  [Ref.  5],  noise  filtering  by 
use  of  local  statistics  [Ref.  6],  and  median  filtering  are  used,  as 
post- processors  to  the  contrast  enhancement  operation.  The  results 
are  shown  in  next  sections.  Another  technique  for  contrast 
manipulation, contrast  enhancement  by  use  of  local  statistics  [Ref. 
6]  was  also  used,  as  an  alternative  to  the  adaptive  filtering 
algorithm.  These  results  are  also  given  and  compared  to  those 
previously  obtained. 

A.     PROCESS  TYPE  1 

This  section  shows  a  combination  of  the  adaptive  filtering  and 
short  space  spatial  subtraction  algorithms.  In  Figure  2.2  two 
degraded  images  are  shown.  The  images  are  characterized  by  poor 
contrast  due  to  the  turbid  water  viewing  conditions.  These  images 
were    processed    for    contrast     enhancement     using     the     adaptive 
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filtering  technique,  but  the  noisy  background  was  heavily 
accentuated  as  shown  in  Figure  2.4.  Then  the  spectral  subtraction 
algorithm  was  used  to  smooth  the  noise  resulting  in  the  images 
shown  in  Figure  4.1.  Inspection  of  these  images  shows  that  this 
particular  algorithm  was  succesful  in  smoothing  the  background 
noise  present  in  the  contrast-enhanced  images,  and,  at  the  same 
time,  most  of  the  information  is  preserved.  Particularly  note,  in 
Figure  4 . 1  (a) ,  how  details  in  the  beam-shaped  object  in  the  center 
of  the  image  are  brought  back,  as  well  as  those  of  the  "fish"  in 
the  lower  right  corner  of  the  image.  The  same  observation  can  be 
made  about  Figure  4 . 1  (b) ,  in  which  the  noise  reduction  was 
achieved,  at  the  cost  of  some  signal  degradation,   in  this  case. 

This  algorithm  requires  the  power  spectral  density  of  the  noise 
to  be  estimated.  In  this  case,  this  was  accomplished  by  taking  a 
portion  of  the  noisy  background,  without  any  objects  and  using  it 
to  estimate  the  spectrum.  Figure  4.2  shows  a  representation  of  the 
estimated  spectrum  of  the  noise  density  thus  obtained.  In  the 
figure,  the  magnitude  of  the  power  spectrum  is  represented  as  the 
intensity  of  a  2-D  image,  each  of  the  image  dimensions  being  the 
corresponding  spatial  frequency  coordinates. 

For  the  short  space  implementation,    it  is  required  to  process 
the   degraded   image   by   dividing   it   into   subimages   and   processing 
each  subimage  separately.    In  dividing  the  degraded  image,   the  size 
of  the  subimage  must  be  such  that  the  image  luminance  in  it  can 
be  approximated  by  a  stationary  random  process.     In  addition,   the 
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(a) 


Figure  4.1. 


(b) 

Images  in  Figure  2.4  Processed  with  the  Short 
Space  Spatial  Subtraction  Algorithm 
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Figure  4.2.     Power  Spectrum  of  the  Noise  in  Figure  2 . 4  (a) 

window  function  nnust  be  a  smooth  function  in  order  to  avoid 
possible  discontinuities  that  may  appear  at  the  subimage  boundaries 
in  the  processed  image.  Thus,  two  types  of  smooth  functions  were 
used-  separable  2-D  triangular  and  Manning  windows-of  size  32x32 
pixels,  overlapped  with  its  neighboring  window  by  half  the  window 
duration  in  each  dimension.  It  was  noted  that  there  was  no 
appreciable  difference  in  using  either  window  -triangular  or 
Manning. 

The  parameter  k,   appearing  in  Equation  3.1  is  a  scaling  factor 
that  normalizes  the  power  spectral  densities,   and  was  evaluated  as 


where  w (1^,12)  is  the  window  function  used, 
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B.  PROCESS  TYPE  2 

In  this  section,  a  combination  of  the  adaptive  filtering  for 
Image  enhancement  and  noise  filtering  by  use  of  local  statistics  is 
presented .  The  images  in  Figure  2 . 4  were  processed  with  the  noise 
filtering  algorithm,  using  different  values  for  the  noise  variance. 
Based  in  the  data  used  for  generating  the  noise  power  spectrum 
shown  in  Figure  4.1,  a  noise  variance  of  60  was  estimated,  for  the 
image  in  Figure  2 . 4  (a) .  For  the  image  in  Figure  2 . 4  (b) ,  a  noise 
variance  of  45  was  estimated.  The  algorithm  was  found  to  be  quite 
sensitive  to  the  size  of  the  window  used  to  evaluate  the  local 
statistics,  the  a  posteriori  mean  and  variance.  Window  sizes  of 
8X8,  16X16,  and  24X24  pixels  were  used  without  noticing  any 
appreciable  improvement  in  the  processed  image .  Figure  4 . 3  shows 
the  results  of  processing  the  images  in  Figure  2 . 4  using  a  window 
size  of  32X32  and  the  noise  variances  mentioned  above.  Note  that 
a  relative  reduction  in  the  background  noise  was  achieved,  with  the 
processed  image  resembling  the  effect  of  low-pass  filtering  in  the 
background,  but  without  blurring  the  edges  of  the  objects  in  the 
images.  When  larger  sizes  for  the  window  were  used,  some 
distortion  in  the  images  begun  to  appear. 

C.  PROCEgS  TYPE  ? 

In   this  section,    the   result   of   applying  a    median    filter    to   a 
contrast-enhanced  image,   using  the  adaptive  filtering  algorithm,    is 
presented.     The    image  in   Figure  2.4(a)    was   processed  with    the 
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(b) 

Figure  4  3     Images  in  Figure  2.4  Processed  with  Noise  Filtering 
by  Use  of  Local  Statistics  Algorithm 
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median  filter  available  in  the  Spider  Library,  using  a  window  size 
of  5x5,  and  making  succesive  iterations  of  the  processed  image, 
until  appreciable  image  degradation  start  to  show  up.  The  image 
obtained  after  8  iterations  is  shown  in  Figure  4.4.  It  can  be 
observed  that  some  noise  smoothing  was  obtained.  Even  when  the 
resulting  image  is  not  noise-free,  the  remaining  noise  is  of  much 
lower  spatial  frequency  than  that  of  the  image  before  the 
median-filtering,   and  thus  seems  less  objectionable. 

Detailed  characteristics  of  the  signal  are  preserved  and  some 
details  have  actually  been  enhanced.  This  is  true  in  the  case  of  the 
fish  in  the  lower  right  corner  of  the  image.  A  comparison  with  the 
use  of  3X3  window  for  the  filter  operation  led  to  the  conclusion 
that  the  smaller  window  size  gives  slightly  more  fidelity  in  the 
signal  preservation.  However,  since  convergence  for  the  3X3 
window  IS  much  slower,  this  choice  is  undesirable.  A  reduction  in 
the  dynamic  range  of  the  median-filtered  image  can  also  be 
observed,  which  presents  an  effect  similar  to  that  of  low-pass 
filtering  in  the  background  region  of  the  image.  However,  this 
procedure  does  not  have  the  undesirable  characteristic  tendency  of 
the  low-pass  process  to  blur  the  edges  of  the  signal. 
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Figure  4.4.     Median  Filtering  of  Image  in  Figure  2.4(a) 

D.     PROCESS  TYPE  4 

In  this  section  the  application  of  the  contrast  enhancement  by 
use  of  local  statistics  algorithm  is  presented.  In  Figures  4.5  and 
4.6  four  low-contrast  noise-degraded  images  are  shown.  All  of 
them  are  512  X  512  pixels  in  size  and  each  pixel  is  represented  by 
8  bits.  The  images  in  Figures  4.5  (a),(b)  and  4.6(a)  are 
characterized  by  a  relatively  dark  background  with  a  diffused 
effect  around  the  objects  in  the  pictures  resulting  in  low  contrast. 
The  image  in  Figure  4.6(b)  has  an  almost  uniform  brightness 
throughout  the  picture,  which  makes  it  difficult  to  distinguish 
between  the  fish  and  beam  shapes  from  the  uniform  background. 

Figure    4.7(a)    corresponds    to    the    image    in    Figure    4.5(a) 
processed  with  the  contrast  enhancement  algorithm,    using  a  value 
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of  4  for  the  gain  k  of  Equation  3.10.  (Recall  that  k  is  the  ratio  of 
the  local  standard  deviation  of  the  processed  image  to  the  original 
standard  deviation) .  Note  the  enhancement  in  contrast  achieved, 
which  IS  particularly  evident  in  the  cable  connecting  the  two 
objects  in  the  lower  part  of  the  image. 

The  image  in  Figure  4 . 5  (b)  was  processed  using  a  value  of 
k=6,  with  results  shown  in  Figure  4.7(b).  The  higher  value  was 
used  to  try  to  enhance  the  fish  shape  in  the  lower  right  part  of 
the  image.   This  enhancement  was  achieved  to  some  extent. 

Figure  4.8(a)  shows  the  result  of  processing  the  image  in 
Figure  4.6(a).  In  this  case,  a  value  for  k  of  4  was  again  used, 
producing  an  image  that  is  a  much  crisper  version  of  the 
corresponding  degraded  image.  Also  note  that  details,  such  as  the 
particular  shape  of  the  object  in  the  picture,  are  enhanced.  The 
noise  in  the  background  is  not  greater  than  that  in  the  original, 
which  means  that  the  desired  result  of  contrast  enhancement  was 
achieved  without  a  significant  noise  increase. 

The  results  of  processing  the  image  in  Figure  4 . 6  (b)  are  shown 
in  Figure  4.8(b).  In  this  case  a  value  of  k=5  was  used,  and  the 
contrast  of  the  two  objects  in  the  picture  are  increased  with 
practically  no  increase  in  the  noise. 
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(b) 
Figure  4,5        Original  Degraded  Images 
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(b) 

Figure  4.6.       Original  Degraded  Images 
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(a) 


SnO    SBS   290 


(b) 

Figure  4.7.       Images  in  Figure  4.5  Processed  with  the  Contrast 
Enhancement  by  Use  of  Local  Statistics  Algorithm 
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(b) 

Figure  4.8.       Images  in  Figure  4.6  Processed  with    The  Contrast 
Enhancement  by  Use  of  Local  Statistics  Algorithm 
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V.  CONCLUSIONS 

Restoration  of  images  degraded  by  turbid  water  was  considered 
in  this  thesis,  Previous  work  on  the  subject  included  the  use  of  an 
adaptive  filtering  algorithm  for  enhancing  the  local  contrast.  This 
algorithm  was  applied  to  images  characterized  by  poor  contrast  due 
to  the  turbid  water  viewing  conditions.  The  technique  yielded  good 
contrast  enhancement,   but  tended  to  accentuate  the  noise. 

In  our  work  some  techniques  for  noise  smoothing  were  used  as 
post-processors  to  the  contrast  enhancement  operation. 

The  specific  techniques  used  were  short  space  spectral 
subtraction,  noise  filtering  by  use  of  local  statistics,  and  median 
filtering.  The  spectral  subtraction  algorithm  yields  very  good  noise 
reduction,  but  tends  to  introduce  some  signal  degradation.  The 
filtering  technique  that  uses  local  statistics  -local  mean  and 
variance-  produced  noise  smoothing  in  the  background,  This 
technique  was  found  to  be  very  sensitive  to  the  size  of  the  window 
in  which  the  local  statistics  were  evaluated.  Median  Filtering  did 
not  eliminate  the  noise  completely,  but  the  remaining  noise  was  of 
much  lower  spatial  frequency  than  the  original  noise  and  thus 
seemed  less  objectionable.  Better  noise  filtering  was  achieved  when 
the  degraded  image  was  median-filtered  iteratively  until  some 
signal  degradation  began  to  appear. 
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As  an  alternative  to  the  adaptive  filtering  algorithnn,  an 
algorithm  for  contrast  enhancement  by  use  of  local  statistics  was 
used.  This  particular  technique  yielded  very  good  results  without 
significant  noise  increase. 

All  processes  implemented  in  this  thesis  proved  to  be  quite 
computationally  intensive.  This  was  especially  true  for  the 
techniques  that  require  the  evaluation  of  local  statistics.  The 
ultimate  goal  for  the  problem  addressed  in  this  and  previous  work, 
is  to  perform  the  enhancement  in  real  time.  This  suggests  a 
natural  continuation  for  the  work  done  so  far  on  this  topic,  in 
order  to  implement  the  algorithms  for  enhancement  in  real  time. 
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APPENDIX 
COMPUTER  PROGRAMS 


PROeRAM  SPATSUB 

C 

c 

C  THIS  PROGRAM  IS  THE  IMPLEMENTATION  OF  THE  SHORT  SPACE 

C  SPECTRAL  SUBTRACTION  ALGORITHM 

C 

C 

c 

BYTE  BYTEIM6(5I2.512) 

CHARACTER*  1  COEN 

CHARACTER* 40  INFILENAME.OUTFILENAME 

INTEGER* 4  IMGBUFF(512,512).  INTARRAY(512,512) 


c  - 

•  *  *  INPUT  FILENAMES  :  INPUT  &  OUTPUT  *  •  * 

TYPE  10 

10 

FORMAT  (•  INPUT  FILENAME  =>  '.$) 

ACCEPT  15,  INFILENAME 

15 

FORMAT  (A40) 

C 

ACCEPT  67.  IWY 

67 

FORMAT  (12) 

TYPE  68 

68 

FORMAT  (•  CaUMN  SIZE  Of  WINDOW  =>  't) 

IWX=16 

C 

ACCEPT  69.  IWX 

69 

FORMAT  (12) 

IERR=0  "' 

•••  INPUT  ORIGINAL  IMAGE  ••* 

CALL  INPUTIMG  (INFILENAME.BYTEIMG.  IMGX.  IMGY) 

*  »  *  CHANGE  IMAGE  FROM  BYTE  TYPE  TO  INTEGER  TYPE  »  *  * 

CALL  BYTE_TOJNTEGER  (BYTEIMG.  IM6BUFF.  IMGX.  IMGY.  lERR) 
IFdERRR  EQ.  1)60  TO  100 

•••  PERFORM  THE  SPATIAL  SUBSTRACTION  **• 
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CALL  CORTO  (1M6BUFF.  INTARRAY,  IM6X,  IMGY.  IWY.  lERR) 

IFdERR  EQ.  1  )  60  TO  100 

C  •  *  »  CHANGE  IMAGE  FROM  INTEGER  TYPE  TO  BYTE  TYPE  *  *  * 

CALL  INTEGEPLTO_BYTE  (INTARRAY.  BYTEIMG.  IMGX,  IM6Y.IERR) 
IFdERR  .EG.  1)  GO  TO  100 

C  *  *  •  OUTPUT  THE  PROCESSED  IMAGE  •  •  • 

CALL  0UTPUTIN6  (OUTFILENAME.  BYTEIMG.  IMGX.  IMGY) 

STOP 

CONTINUE 

TYPE   *.■  II  FORCED  TO  EXIT  !!  • 

STOP 

END 

SUBROUTINE  CORTO  (ENTRADA.  SALIDA.m6X.m6Y.IWY.IERR) 

C  THIS  SUBROUTINE  PERFORMS  THE  SPATIAL  SUBTRACTION  ALGORITHN 

C  USING  THE  SHORT  SPACE  IMPLEMENTATION  TECHNIQUE  FOR  IMAGE 

C  RESTORATION 

INTEGER* 4  ENTRADA  (IMGX.IMGY).  SALIDA(IMGX,IMGY).  INTE(32.32) 

REALM    AR(32.32),AI(32.32).BR(32.32),BI(32.32). 
+       ST(64).    CT(64).    LBR(64).  BRR(32.32).  WIND0W(32.32). 
+       R£NTRADA(512.5I2).  IMGWIN(32.32).  IARRAY(32.32).  MAG(32.32). 
+       PHAS(32,32).PSD(32.32) 

DO  I   -  I.IMGX 
DO  J  -  I  .IMGY 

RENTRADA  (I, J)  =  FLOATJ(ENTRADA(I.J)) 
END  DO 
END  DO 

TYPE  * /INTEGER  IMAGE  .CONVERTED  TO  REAL  DATA' 
C  CALL  MCHECK(RENTRADA.IMGX) 

C  •••   BUILD  UP  THE  WINDOW  FUNCTION  ••  • 

CALL  BLTWINDOW  (16,WIND0V/) 
C  TYPE  *  .WINDOW  FUNCTION  BUILF 

C  CALL  MCHECK(WINDOW.32) 
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AK  -0.0 

DO  1  =  1 ,32 
DO  J  -1,32 

AK  =  AK+(WIND0W(IJ)««2 
END  DO 
END  DO 

TYPE  *,AK->...  -.AK 

C  *  •  •  INITIALIZE  OUTPUT  ARRAY  •  •  * 

DO  L  =  1  .IM6 
DO  M  -  1  ,IM6Y 

SALIDA  (L.M)  =  0 
END  DO 
END  DO 

C  TYPE  ".OUTPUT  ARRAY  INITIALIZED' 

C  CALL  MCHECK(SALIDA.iri6X) 

C  *  *  *  ESTIMATE  POWER  SPECTRUM  DENSITY  OF  THE  NOISE  *  *  « 

CALL  ESTINOISE(RENTRADA.PSD) 

C  TYPE  •  ,POWER  SPECTRUM  DENSITY  OF  THE  NOISE  ESTIMATED' 

C  CALL  MCHECK(PSD.32) 

C  •  •  *  MAIN  PROCESS  SECTION  *  •  • 

C  INSERT  LIMITS  HERE 

DO  1-1.31 
DO  J=  1.31 

LL  -((l-1)*16)+1 
MM  -  ((J-1)*16)+1 
DO  L- 1.32 
DOM-  1,32 

IA-LL+L-1 
B  -MM+M-1 
C  IF(IA  .6E.  16.ANDIA  ,LE.  112  .AND. 

C  +  IB  .6E.  16  .AND.  IB  .LE.  112)  THEN 

IM6WIN(L.M)-RENTRADA(IA.IB) 
IARRAY(L,M)-WIND0W(L.M)»IM6WIN(L.M) 

C  ELSE  IF(I  A  .6E .  1 6  .AND .  I A  .LE .  1 1 2)  THEN 
C  lARRAY(L.M)-  0.5»R£NTRADA(IA.IB) 

e  ELSE  IF(IB  .6E.  16  AND,  IB  LE.  1 12)  THEN 
C  IARRAY(L,M)-0.5»RENTRADA(IA.IB) 
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C  ELSE 

C  lARRAY  (L.M)=RENTRADA(IA.IB) 

C 

C  END  IF 

END  DO 

END  DO 

C  TYPE  *  ."DATA  ARRAY  HAS  BEEN  WINDOWED' 

C  CALL  MCHECK  (IARRAY.32) 

C  *  "  *  TRANSFORM  DATA  ARRAY  INTO  COMPLEX  FORMAT  *  *  * 

DO  L=  1.32 
DO  M=  1.32 

AI(L.M)  =  0. 

AR(L.M)  =  lARRAY(L.M) 
END  DO 
END  DO 

C  TYPE *. CHECK  REAL  PART' 

C  CALL  MCHECK  (AR.32) 

C  TYPE  ». CHECK  IMAGINARY  PART  (ZEROS)' 

C  CALL  MCHECK  (AI.32) 

C  ••  *  *  EVALUATE  2-0  FFT  OF  DATA  IMAGE  "  *  » 

CALL    FFTS2(AR.AI.BR.BI.32.32.ST.CT.LBR.64.2.JERR) 
C  CALL  MCHECK  (BR.32) 

DO  IX  =  1.32 
DOIY=  1.32 

BRRCIX.IYKABS    (CMPLX(BR(IX.IY).BI(IX.IY))) 
END  DO 
END  DO 

C  TYPE  •.'FFT  DONE. CHECK  MAGNITUDE' 

C  CALL  MCHECK  (BRR.32) 

C  CALL  DISPECT  (BRR) 

C  »  •  •  PERFORFM  SPATIAL  SUBSTRACTION  »  •  • 

ALFA  =  1 .5 
C  PSD  =  40.0 

DO  IX  =  1 .32 
DO  lY  =  1 ,32 

BR(IX.IY)  =  ((BRR(IX.IY))»»2)    -    ((ALFA/AK)  *  ((PSD(IX.IY)  »»2))) 
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IF  (BRdX.IY)  IT.  .0.)  THEN 

MAG  (IX. lY)  =  0 
ELSE 

MAGdX.IY)  =  SQRKBRdX.IY)) 
END  IF 
END  DO 
END  DO 

C  TYPE    -  .-SPATIAL  SUBT  DONE  CHECK  RESULT' 

C  CALL  MCHECK  (MAG) 

C  •  *  •  EVAL ,  MAGNITUDE  OF  INV.  FFT  •  *  • 

DOIX=  1.32 

DOIY=  1,32 

BR(IX.IY)  =  MAG(IX.IY)«COS(PHAS(IX.IY) 
Bl  (IX. lY)  =  MAG(IX.IY)*SIN(PHAS(IX.IY) 

END  DO 

END  DO 

CALL    FFTS2(BR.BI.AR.AI.32.32.ST.CT,LBR.64,-2.JERR) 

C  TYPE", 'INV.  FFT  DONE.  CHECK  REAL  PART' 

C  CALL  MCHECK(AR) 

C  TYPE  ".CHECK  IMAGINARY  PART  (ZEROS)  ' 

C  CALL  MCHECK(AI) 

C  •  »  "  BACK  TO  INTEGER  FORMAT  »  •  » 

DO  IX  =  1 .32 
DOIY=  1.32 

INTE(IX.IY)  =  JINT(AR(IX.IY)) 
END  DO 
END  DO 

C  TYPE", INTEGER  CONVERSION  DONE.CHECK  IT' 

C  CALL  MCHECK  (INTE) 

C  •  •  *  BUILD  UP  THE  OUTPUT  ARRAY  •  •  • 

DOIX=  1.32 
DO  lY  =  1 .32 

LL  =((I-1)*16K1 

MM((J-1)*16>+1 

IXX=LL+IX-1 

lYY  =MM+IY-1 

SALIDA  (IXX.IYY)  =  SALIDA(IXX.IYY>+INTE(IX.IY) 
C  TYPE".1XX  AND  lYY  =>',IXX.IYY 
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END  DO 
END  DO 

END  DO  !!  END  DO  J 
END  DO  !!  END  DO  I 

100  TYPE».T1A1N  PROCESS  DONE' 

C  CALL  TRUNCATE(5ALIDA,IM6S.IM6Y) 

CALL   SCALE^DJUST(SALIDA.IM6X.IMGY.0.0.255,0) 
TYPE    *. SPATIAL  PROCESS  DONE' 

RETURN 

END 


SUBROUTINE  BLTWIND0W(IWY.NiN2) 

REALM  Nl(32),    NIN2(32,32).  N2(32,32) 


L=0 

PI  o4-ATAN(l.) 

DOI+  1.32 

II  =1-1 

N1(I)  =  0.5-(0,5*COS(2*PI»II/31)) 

N2(I)  =  N1(I) 

c 

TYPE  •.NKI).  N2(l) 

c 

IF(I  .6T.  16)  THEN 

c 

L=L+1 

c 

N1(l)  =  l  -(2«L) 

c 

ELSE 

c 

N1(l)=l 

c 

END  IF 

END  DO 

L=0 

DO  1  =  1  .  32 

DOJ=  1.32 

c 

IF  (J  .6T.  16)  THEN 

c 

L=L+1 

c 

NIN2(U)  -  NKI)  *  (J-2'»1) 

c 

ELSE 

c 

NIN2(I,J)  =  NKI)  »  J 

c 

END  IF 

NIN2(U)  '  N1(I)*N2(J) 

c 

TYPE  •.  NIN2(I.J) 

END  DO 

L=0 
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END  DO 
RETURN 
END 

SUBROUTINE  ESTINOISE  (RENTRADA.PSD) 


C  THIS  ROUTINE  ESTIMATES  THE  POWER  SPECTRAL  DENSITY  Of 

C  THE  NOISE  IN  AN  IMA6E.  BY  TAKING  THE  MAGNITUDE  (SQUARED) 

C  Of  THE  DISCRETE  SPACE  FOURIER  TRANSFORM  Of  A  (NOISY)  PART 

C  Of  THE  BACKGROUND. 

C  TO  CHOOSE  APPROPIATE  REGION.  SET  PARAMETERS  I  AND  J 

C  ACCORDINGLY. 

REALM  RENTRADA(512,512).  PSD(32.32).  AR(32.32).  Al(32.32). 
BR(32.32).  Bl(32.32).  ST(32.32).  CT(32.32).  LBR(64) 

C  •  •  "  TAKE  A  PORTION  Of  THE  IMAGE  BACKGROUND  •  •  • 

J-80 
K-80 

DO  L- 1.32 
DO  M=  1,32 

LL-L+J 

MM=M+K 

AR(L.M)-RENTRADA(LL.MM) 
END  DO 
END  DO 

C  •  •  •  INITIALIZE  IMAGINARY  PART  Of  DATA  •  •  » 

DO  L- 1,32 
DO  M=  1,32 

AI(L.M)-0. 
END  DO 
END  DO 

C  "  •  •  EVALUATE  MAGNITUDE  Of  2-0  f FT  "  •  • 

CALL    FFTS2(AR.AI.BR.BI.32.32.ST.CT.LBR.64.2.JERR) 

DO  L- 1.32 
DOM- 1.32 

PSD(L.MKABS(CMPLX(BR(L.M).BI(L.M))) 
END  DO 
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END  DO 
RETURN 
END 

SUBROUTINE    INPUTin6(INFILENAnE.BYTEin6.in6X.m6Y) 

C  READ  IN  INPUT  IMAGE  DATA  IN  BYTE  TYPE 

C  ASSUME  INPUT  DATA  FILE  IS  A  DIRECT  ACCESS  FILE 

BYTE    BYTEIM6(IM6Y,IM6X) 

CHARACTER* 40  INFILENAME 


0PEN(UNIT=2,  FILE=INFILENAME.  STATUS=aD'.  ACCESS=DIR£CT' 
+       RECORDTYPE=TIXED') 


DO  I=1.IM6Y 

READ(2'I)  (BYTEIMGd.  J),  J=1,  IM6X) 

END  DO 

TYPE  *.  •  *  INPUT  IMAGE  HAS  BEEN  READ  " 

CL0SE(2) 

RETURN 
END 


SUBROUTINE  BYTE_T0_INTE6ER  (BYTEDATA.  INTDATA.  IX.  lY.  lERR) 

C 

C  THIS  PROGRAM  CHANGE  DATA  (USUALLY  2D  IMAGE  DATA)  IN  THE  BYTE  DATA 

C  TYPE  INTO  THE  INTEGER  TYPE 

C 

BYTE       BYTEDATA  (lY.  IX) 
INTE6ER-4  INTDATA  (lY.  IX) 

DO  I  =  1 .  lY 
DOJ=  1.IX 

IF  (BYTEDATA  (I. J)  .GE.  -128 
+        AND.  BYTEDATA  (U)  IT.  0)  THEN 

INTDATAdJ)  =  BYTEDATA  (I. J)  +  256 

ELSE  IF  (BYTEDATA(U)  GE.O 
+       .AND.  BYTEDATAdJ)  LE.  127)  THEN 

INTDATA(U)  =  BYTEDATAdJ) 

ELSE 
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TYPE*;?JDATA  OUT  OF  RANGE!  ROW,  Ca.  DATA  =>  •, 
+       I.  J.  BYTEDATA(I.J) 
IERR  =  1 
GOTO  100 
ENDIF 

END  DO 
END  DO 

TYPE*.  •  BYTE  TO  INTEGER  CONVERSION  HAS  BEEN  DONE' 

100  CONTINUE 

RETURN 
END 

SUBROUTINE  INTE6ER_T0_BYTE(INTDATA.  BYTEDATA.  IX.  lY.  lERR) 

C 

C  THIS  PROGRAM  CHANGE  DATA  (USUALLY  2D  IMAGE  DATA)  IN  THE  BYTE  DATA 

C  TYPE  INTO  THE  INTEGER  TYPE 

C 

BYTE  BYTEDATAdX.  lY) 
INTEGER*  4  INTDATAdX.  lY) 

DO  1=1, IX 
DO  J=1.IY 

IF  (INTDATAdJ)  IT.  0  .OR.    INTDATAdJ)  .6T.  255)  THEN  TYPE*.'?? 
DATA  OUTOF  RANGE!  ROW.Ca.DATA  =>  '. 
+       I.  J.  INTDATAdJ) 
lERR  =  1 
GOTO  100 

ELSE  IF  (INTDATAdJ)  .IE.  127)  THEN 
BYTEDATAdJ)  =  INTDATA(IJ) 
ELSE 

BYTEDATAdJ)  =  INTDATAd J)  -  256 
ENDIF 

END  DO 
END  DO 

TYPE  • ;  INTEGER  TO  BYTE  CONVERSION  HAS  BEEN  DONE' 

100  CONTINUE 
RETURN 
END 
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PR06RAn  LOC_STAT 


C  THIS  PROGRAM  IS  THE  IMPLEMENTATION  OF  THE  IMAGE  ENHANCEMENT  AND 

C  NOISE  FILTERING  ALGORITHM 

C 

C 

BYTE  BYTEIMG  (512.512) 

CHARACTER*  1  COEN 

CHARACTER«40  INFILENAME.  OUTFILENAME 

INTEGER* 2  IM6BUFF(512.512),  LAVEIM6(512.512).  L0CMEAN(512.512). 
-       INTARRAY(512.512) 


C 

**•  INPUT  FILENAMES  :  INPUT  &  OUTPUT  •** 

TYPE  10 

10 

FORMAT  (•  INPUT  FILENAME  ->  '.$) 

ACCEPT  15.  INFILENAME 

15 

FORMAT  (A40) 

TYPE  20 

20 

FORMAT  (■  OUTPUT  FILENAME  »  \%) 

ACCEPT  25.  OUTFILENAME 

25 

FORMAT  (A40) 

C 

•"*  INPUT  ORIGINAL  IMAGE  SIZE  •»* 

TYPE  50 

50 

FORMAT  (•  ROW  SIZE  OF  IMAGE  ->  '.$) 

ACCEPT  55.  IMGY 

55 

FORMAT  (14) 

TYPE  60 

60 

FORMAT  (•  CaUMN  SIZE  OF  IMAGE  =>  *.$) 

ACCEPT  65,  IM6X 

65 

FORMAT  914) 

C 

•»*  INPUT  WINDOW  SIZE  *•* 

TYPE  66 

66 

FORMAT  (■  ROW  SIZE  OF  WINDOW  =>  •.$) 

ACCEPT  67,  IWY 

67 

FORMAT  (12) 

TYPE  68 

68 

FORMAT  (•  CaUMN  SIZE  OF  WINDOW  =>  ',$) 

ACCEPT  69,  IWX 

69 

FORMAT  (12) 

C 

*  *  *  SELECT  CONTRAST  ENHANCEMENT  OR  NOISI 

TYPE  70 
70  FORMAT  (  'CONTRAST/NOISEFIL  ?  (C/N)  =•>•.$) 
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ACCEPT  75.  CE 
75  FORMAT  (AD 


IERR  =  0 
C  *  *  •  INPUT  ORIGINAL  IMAGE  •  •  • 

CALL  INPUTIM6  (INFILENAME.  BYTEIM6.  IM6X.  IMGY) 

C  *  •  •  CHANGE  IMAGE  FROM  BYTE  TYPE  INTEGER  TYPE  •  •  « 

CALL  BYTE.TOJNTEGER  (  BYTEIMG.  IM6BUFF.  IMGX.  IMGY,  lERR) 
IFdERR  EQ.DGOTO  100 

C  •  •  •  CALCULATE  LOCAL  MEAN  ARRAY  •  •  • 

CALL  LOCALJIEANdMGBUFF.  LOCMEAN,  IMGX.  IMGY.  IWX.  IWY.  lERR) 
IFdERR  EQ.I  )60T0  100 

IF  (CE  EQ.  -C  .OR.  CE  .EQ.  c)  THEN 

C  •  •  •  INPUT  K  =  STANDARD  DEVIATION  RATIO  •  "  * 

TYPE  80 
80  FORMAT  C  STAND.  DEV.  RATIO  K=>'.  $) 

ACCEPT  85.  K 
85  FORMAT  (14) 

CALL  ENHACONT  (IM6BUFF.  LOCMEAN.  INTARRAY.  IMGX.  IMGY.  K.  lERR) 

ELSE 

C  •  ••INPUTNV  =  NOISE  VARIANCE  •»• 

TYPE  90 
90  FORMAT  C  NOISE  VARIANCE  NV  => ".  t) 

ACCEPT  95 
95  FORMAT  (M) 

CALL  NOISE_f  ILT  (  IM6BUFF.  LOCMEAN.  INTARRAY.  NV.  IMGS.  IMGY.  lERR) 
ENDIF 

SCLMIN  -  0.0 
SCLMAX  =  255.0 
C  CALL  SCALE_ADJUST  (INTARRAY.  IMGX.  IMGY.  SCLMAX.  SCLMIN) 

CALL  TRUNCATE  (INTARRAY.IMGX.IMGY) 

C  *  •  *  CHANGE  IMAGE  FROM  INTEGER  TYPE  TO  BYTE  TYPE  •  «• » 

CALL  INTE6ER_T0_BYTE  (  INTARRAY.  BYTEIMG.  IMGX.  IMGY.  lERR  ) 
IFdERR  EQ.  DGOTO  100 

C  *  •  •  OUTPUT  THE  PROCESSED  IMAGE  •  •  - 
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CALL  OUTPUTine  (  OUTFILENAME.  BYTEIMG,  IMOX.  IM6Y  ) 


STOP 
100  CONTINUE 


TYPE  V  !!  FORCED  TO  EXIT  I! 
STOP 

END 


SUBROUTINE   ENHACONT   (in6BUFF.L0CnEAN.INTARRAY.in6X.in6Y.IC.IERR) 


C  THIS  SUBROUTINE  IMPLEMENTS  THE  CONTRAST  ENHANCEMENT  BY  USE 

C  OF  LOCAL  STATISTICS  ALGORITHM. 

C 


INTEGER* 2  IMGBUFFCIMGY.IMGX)  .  LOCMEANdMGY.IMGX)  . 
INTARRAYdMGY.IMGX) 

»  •  •  EVALUATE  LOCAL  CONTRAST  AND  MODIFY  IT  WITH  K  »  *  * 

DO  I  -  1  IM6Y 

DO  J  -  I.IMGX 

INTARRAY(IJ)  -  K«(IM6BUFF(I,J)  -  LOCMEANdJ)) 

END  DO 

END  DO 

•  -  •  ADD  MODIFIED  LOC.  CONTRAST  AND  LOCAL  MEAN  •  •  • 

DO  I  -  1  .IM6Y 

DOJ  =  I.IMGX 

INTARRAY  (I, J)  =  INTARRAY  (I, J)  +  LOCMEAN  (I.J) 

END  DO 

END  DO 

•  *  *  CORRECT  ANY  OVERFLOW  »  »  * 

DO  I  -  1  .IM6Y 

DO  J  =  1  .IMGX 

INTARRAY(U)  =  K*(IMGBUFF(I.J)  -  LOCMEAN(I.J)  ) 

END  DO 

END  DO 

•  •  •  ADD  MODIFIED  LOC.  CONTRAST  AND  LOCAL  MEAN  *  *  • 

DO  I  =  1  .IMGY 

DO  J  =  I.IMGX 

INTARRAY  (I.J)  -  INTARRAY  (I.J)  +  LOCMEAN  (l,J) 

END  DO 
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END  DO 

»  "  •  CORRECT  ANY  OVERFLOW  •  •  • 

SCLMIN  =  0.0 

SCLMAX  =  255,0 

CALL  SCALE_ADJUST  (INTARRAY.  IMOX.  IM6Y.  SCLMAX.  SCLMIN) 

RETURN 
END 


SUBROUTINE    N0iSE-FILT(m6BUFF.L0CnEAN.INTARRAY.NV.m6X. 
*       IMGY.IERR) 


C  THIS  SUBROUTINE  IS  THE  IMPLEMENTATION  OF  THE  NOISE  FILTERING 

C  ALGORITHM.  CALLED  BY  THE  NOISE  FILTERING  AND  CONTRAST  ENHANCEMENT 

C  BY  USE  OF  LOCAL  STATISTICS  PROGRAM. 

C 

C  REFERENCES; 

C 

c 

INTEGER* 2   IMGBUFFdMGY.IMGX).  LOCMEAN(IMGY.IMGX). 
+       INTARRAYdMGY.IMGX) 

REAL-8  AK(128.128)  .0(128.128) 

C  *  *  *  EVALUATE  THE  VARIANCE  OF  THE  IMPUT  IMAGE  »  *  * 

TYPE  -.CHECK  IMGBUFF  ARRAY' 
C  CALL  MCHECK(IMGBUFF) 

TYPE  •  .CHECK  LOCMEAN  ARRAY  ' 
C  CALL  MCHECK(LOCMEAN) 

DO  I  =  1  .IMGY 

DO  J  =  1.IMGX 

IDIF  =  IMGBUFF(I.J)  -  LOCMEAN(I.J) 
C  TYPE  16.  IDIF 

16  FORMAT  (■  DIF  =>■.  MO) 
Q(I.J)  =  (IDIF»*2)-NV 

C         TYPE  17.  Q(U) 

17  FORMATC  Q  =>  •.F8.2) 
END  DO 

END  DO 

TYPE  *.' CHECK  THE  Q  ARRAY  ' 
C  CALL  MRCHECK(Q) 

C  »"*  EVALUATE  THE  GAIN  AK  **  • 

DO  1=  1.IMGY 
DO  J  =  1 . IMGX 
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DEN  =  GKI.J)  +  NV 

IF  (DEN  .EQ.  0.0)  THEN 
DEN  =  1 .0 

ENDIF 

AK(I.J)  =  (GKIJ)  /  DEN) 
END  DO 
END  DO 

**•  ESTIMATE  THE  RESTORED  IMA6E  *** 

DO  I  =  1 . IM6Y 

DO  J  =  1 .  IM6X 

DIP   =  AK(I.J)»{IMGBUFF(I.J)  -LOCMEAN(I.J)) 
INTARRAY(I.J)  =  INT(LOCMEAN(I.J)  +  DIF) 

END  DO 

END  DO 

RETURN 
END 


SUBROUTINE  SCALE-ADJUST  (IMGBUFF.  mCX.  mCY.  SCLMAX.  SCmiN) 

INTEGER* 4  IMGBUFF  (IM6X.  IMGY) 

MAXVAL  =  IMGBUFF  (1.1) 
MINVAL  =  IMGBUFF  (1,1) 

DO  1=1.  IMGY 
DO  J=  1.  IMGX 

IF  (IMGBUFF(I.J)  GT.  MAXVAL)  THEN 
MAXVAL  =  IMGBUFFdJ) 

ENDIF 

IF  (IMGBUFFdJ)  .IT.  MINVAL)  THEN 

MINVAL  =  IMGBUFF(I.J) 
ENDIF 

END  DO 
END  DO 

C  TYPE  10.  MAXVAL.  MINVAL 

CIO  FORMATCMAX  &MIN  VALUES  =>'  218) 

TYPE  ".riAXVAL^   .MAXVAL.  T1INVAAL=   .MINVAL 

INTVAL  =  MAXVAL  -  MINVAL 

SLOPE  =  (SCLMAX  -  SCLMIN)/REAL  (INTVAL) 
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DO  I  =  1 ,  II16Y 
DO  J  =  1 , IMGX 


IF  (IM6BUFF(I.J)  .EQ.  MINVAL)  THEN 
II16BUFF  (I.J)  -  SCLMIN 

ELSE  IF  (IM6BUFF(I.J)  .EQ.  MAXVAL)  THEN 

IMGBUFF  (I.J)  -  SCLMAX 
ELSE 


IDIF  =  II16BUFF  (I.J)  -  MINVAL 

INGBUFFdJ)  =  INT(S)  OPE" IDIF)  +  INT(SCLMIN 


ENDIF 


END  DO 
END  DO 

RETURN 
END 


SUBROUTINE  TRUNCATE(m6BUFF.  IMGX.  IM6Y) 

TRUNCATE  IMAGE  GRAY  VALUES  INTO  0-255 
INTEGER"  4  IMGBUFFdMGY,  IMGX) 

DO  I  -  1 . IMGY 

DO  J-  1.IMGX 

IF  (IMGBUFF(U)  .LT.  0)  IMGBUFF(IJ)  =  0 

IF  (IMGBUFFdJ)  .GT.  255)  IMGBUFF(l.J)  -  255 

END  DO 

END  DO 

RETURN 
END 
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PR06RAn  HEDIFILT 


C  THIS  PROGRAM  IS  THE  IMPLEMENTATION  Of  THE  MEDIAN 

C  FILTERING  IN  NOISE  REDUCTION 


BYTE  BYTEIM6(512.  512) 

CHARACTER* 40  INFILENAME.  OUTFILENAME 

INTEGER* 4  IMGBUFF(512.512).  IHST(256).  INTARRAY(512.  512) 


c 

••*  INPUT  FILENAMES  ;  INPUT  &  OUTPUT  ••• 

TYPE  10 

10 

FORMAT  CINPUT  FILENAME  =>  \t) 

ACCEPT  15.  INFILENAME 

FORMAT  (A40) 

TYPE  20 

20 

FORMAT  ('OUTPUT  FILENAME    =>  '.$) 

ACCEPT  25.  OUTFILENAME 

25 

FORMAT  (A40) 

C 

***  INPUT  ORIGINAL  IMAGE  SIZE  "•* 

C 

TYPE  50 

50 

FORMAT  CROW  SIZE  OF  IMAGE  =>  '.$) 

IMGY=512 

C 

ACCEPT  55.  IMGY 

55 

FORMAT  (14) 

C 

TYPE  60 

60 

FORMAT  (•  CaUMN  SIZE  OF  IMAGE  =>',$) 

IMGX  =  512 

C 

ACCEPT  65.  IMGX 

65 

FORMAT  (14) 

TYPE  *. IMGX. IMGY 

IERR  =  0 

»*«  INPUT  ORIGINAL  IMAGE  ••* 

CALL  IMPUTIMG  (INFILENAME.BYTEIMG.  IMGX,  IMGY) 

«••  CHANGE  IMAGE  FROM  BYTE  TYPE  TO  INTEGER  TYPE  *-• 

CALL  BYTE_TOJNTEGER  (BYTEIMG.  IMGBUFF.  IMGX.  IMGY.  lERR) 
IFdERR  .EQ.  1  )  GO  TO  100 
*  »  *   PERFORM  THE  MEDIAN  FILTERING  •  •  • 
DO  I  =  1 .  20 

CALL  MEDI  (IMGBUFF.  INTARRAY.  IMGX.  IMGY.  5.  5.  IHST.  256) 
DOIX=  1.  IMGX 
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DOIY=  1,  IMGY 

IMGBUFFdX.IY)  =  INTARRAYCIX.  lY) 
END  DO 
END  DO 
END  DO 
C  •  "  •  CHANGE  IMAGE  FROM  INTEGER  TYPE  TO  BYTE  TYPE  •  *  • 

C  TYPE  •.  IMGX.  IMGY 

CALL  INTEGER_TO_BYTE  (INTARRAY.  BYTEIMG.  IMGX.  IMGY,  lERR) 
IFdERR  .EQ.  1)  GO  TO  100 

C  »  *  *  OUTPUT  THE  PROCESSED  IMAGE  *  "  • 

CALL  OUTPUTIMG  (OUTFILENAME.  BYTEIMG.  IMGX.  IMGY) 

STOP 

100       CONTINUE 

TYPE  •.'  !  !  FORCED  TO  EXIT  !  C 
STOP 

END 
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